knitr::opts_chunk$set(fig.path = 'figs/',
echo = TRUE, message = FALSE, warning = FALSE)
library(sf)
library(raster)
library(data.table)
source('https://raw.githubusercontent.com/oharac/src/master/R/common.R')
dir_git <- here()
source(file.path(dir_git, '_setup/common_fxns.R'))Examine biodiversity risk by taxonomic group.
spp_maps <- read_csv(file.path(dir_data,
sprintf('spp_marine_maps_%s.csv', api_version)),
col_types = 'ddcicccl') %>%
filter(comp_assessed == TRUE) %>%
mutate(dbf_file = basename(dbf_file)) %>%
select(iucn_sid, dbf_file)
spp_info <- read_csv(file.path(dir_o_anx,
sprintf('iucn/spp_info_from_api_%s.csv', api_version))) %>%
select(-sciname)
taxa_list <- spp_maps %>%
left_join(spp_info, by = 'iucn_sid') %>%
select(iucn_sid, kingdom, phylum, class, order, family, dbf_file)
fam_list <- taxa_list %>%
group_by(phylum, class, order, family, dbf_file) %>%
summarize(n_spp = n())
order_list <- taxa_list %>%
group_by(phylum, class, order, dbf_file) %>%
summarize(n_spp = n())
class_list <- taxa_list %>%
group_by(phylum, class, dbf_file) %>%
summarize(n_spp = n())shp_to_taxa <- read_csv(file.path(dir_setup, 'raw', 'shps_to_taxa.csv'))
taxa_sums <- data.frame(
unweighted = list.files(file.path(dir_o_anx, 'taxa_sums_mol_2019'),
pattern = sprintf('cell_sum_comp_%s.csv', api_version),
full.names = TRUE),
rr_weighted = list.files(file.path(dir_o_anx, 'taxa_sums_mol_2019'),
pattern = sprintf('cell_sum_comp_rrweight_%s.csv',
api_version),
full.names = TRUE)) %>%
mutate(shp_base = basename(unweighted),
shp_base = str_replace_all(shp_base, '_cell_sum.+|_part.+', '')) %>%
left_join(shp_to_taxa, by = 'shp_base')
# write_csv(taxa_sums %>% select(shp_base, taxon) %>% distinct(),
# file.path(dir_setup, 'raw', 'shps_to_taxa_raw.csv'))
spp_maps <- read_csv(file.path(dir_data, sprintf('spp_marine_maps_%s.csv', api_version)),
col_types = 'ddcicccl') %>%
filter(comp_assessed == TRUE) %>%
mutate(shp_base = tolower(dbf_file) %>%
basename() %>%
str_replace_all(., '_part.+|.dbf', '')) %>%
select(iucn_sid, shp_base) %>%
distinct() %>%
group_by(shp_base) %>%
summarize(n_spp = n()) %>%
left_join(taxa_sums, by = 'shp_base') %>%
select(taxon, n_spp) %>%
distinct() %>%
group_by(taxon) %>%
summarize(n_spp = sum(n_spp)) %>%
arrange(desc(n_spp))reload <- TRUE ### reload applies to building rasters, not the plot
rast_base <- raster::raster(file.path(dir_spatial, 'cell_id_mol.tif'))
land_poly <- sf::read_sf(file.path(dir_spatial, 'ne_10m_land',
'ne_10m_land_no_casp.shp')) %>%
st_transform(crs(rast_base))
taxa_gps <- c("Cnidaria",
"Tracheophyta",
"Chordata: Actinopterygii",
"Chordata: Mammalia",
"Chordata: Chondrichthyes",
"Chordata: Aves",
"Chordata: Reptilia",
"Echinodermata")
map_list <- vector('list', length = length(taxa_gps)) %>%
setNames(taxa_gps)
for(i in seq_along(taxa_gps)) { ### i <- 6
taxon_gp <- taxa_gps[i]
taxon_txt <- taxon_gp %>%
tolower() %>%
str_replace_all('chordata: |[^a-z]+', '')
taxa_info <- taxa_sums %>%
filter(taxon == taxon_gp)
taxa_rast_file <- file.path(dir_o_anx, 'taxa_rasts_mol_2019',
paste0(tolower(taxon_txt), '_rr_comp.tif'))
if(!file.exists(taxa_rast_file) | reload == TRUE) {
cat_msg('Building: ', taxa_rast_file)
rr_weighted_files <- taxa_info$rr_weighted %>% unique()
taxon_df_rr <- parallel::mclapply(rr_weighted_files, mc.cores = 12,
FUN = function (x) { ### x <- rr_weighted_files[1]
read_csv(x, col_types = 'ddd___d___')
}) %>%
bind_rows() %>%
group_by(cell_id) %>%
summarize(mean_risk_sum = 1/sum(sr_rr_risk) * sum(mean_risk * sr_rr_risk),
pct_threat_sum = sum(sr_rr_threatened) / sum(sr_rr_risk)) %>%
select(cell_id, mean_risk_sum, pct_threat_sum)
mean_rr_rast <- raster::subs(x = rast_base, y = taxon_df_rr,
by = 'cell_id', which = 'mean_risk_sum')
raster::writeRaster(mean_rr_rast, filename = taxa_rast_file, overwrite = TRUE)
} else {
cat_msg('File exists: ', taxa_rast_file)
mean_rr_rast <- raster(taxa_rast_file)
}
mean_rr_df <- mean_rr_rast %>%
aggregate(fact = 2) %>%
### should look OK for tiny maps?
rasterToPoints() %>%
as.data.frame() %>%
setNames(c('long', 'lat', 'value'))
cat_msg('Building plot for ', taxon_gp)
### using geom_tile so I can add a size to pixels, thought it is
### significantly slower...
x <- ggplot(mean_rr_df) +
geom_tile(aes(long, lat, fill = value, color = value),
size = .5,
show.legend = FALSE) +
geom_sf(data = land_poly, aes(geometry = geometry),
fill = 'grey96', color = 'grey40', size = .10) +
ggtheme_map(base_size = 9) +
theme(plot.margin = unit(c(.05, .10, .05, .05), units = 'cm'),
panel.background = element_rect(fill = 'grey80', color = NA),
axis.title = element_text(face = 'bold', hjust = 1,
vjust = 0, size = 7,
margin = margin(t = 0, r = 0,
b = 0, l = .1, unit = 'cm')),
axis.title.x = element_blank()) +
coord_sf(datum = NA) + ### ditch graticules
scale_fill_gradientn(colors = risk_cols,
values = risk_vals,
limits = c(0, 1),
labels = risk_lbls,
breaks = risk_brks) +
scale_color_gradientn(colors = risk_cols,
values = risk_vals,
limits = c(0, 1),
labels = risk_lbls,
breaks = risk_brks) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
labs(y = taxon_gp)
map_list[[i]] <- x
} ### set up a dataframe of values to craft a color bar using geom_segment
colorbar_df <- data.frame(x = seq(0, 1, .001), y = -1)
colorbar <- ggplot(colorbar_df) +
ggtheme_plot(base_size = 7) +
geom_segment(aes(x = x, xend = x, color = x),
y = 0, yend = 1, size = 1,
show.legend = FALSE) +
scale_color_gradientn(colors = risk_cols, values = risk_vals,
limits = c(0, 1),
labels = risk_lbls, breaks = risk_brks) +
theme(axis.text.y = element_blank(),
axis.title = element_blank(),
panel.grid.major.y = element_blank(),
plot.margin = unit(c(.05, 0, .05, 0), units = 'cm')) +
scale_x_continuous(labels = risk_lbls, breaks = risk_brks, limits = c(0, 1),
expand = expand_scale(.25)) +
scale_y_continuous(limits = c(0, 1))cat_msg('done, putting together the cowplot')
y <- cowplot::plot_grid(plotlist = map_list, align = 'hv', ncol = 2)
z <- cowplot::plot_grid(y, colorbar,
rel_heights = c(20, 1),
align = 'v', ncol = 1)
plot_w <- 4.75
plot_h <- 6.1
taxa_map_file <- file.path(dir_git, 'ms_figures/fig2_taxa_maps_mol.tiff')
taxa_map_file <- file.path(dir_git, 'ms_figures/fig2_taxa_maps_mol.png')
ggsave(plot = z,
filename = taxa_map_file,
width = plot_w, height = plot_h, units = 'in', dpi = 300)spp_ranges <- read_csv(file.path(dir_data,
sprintf('iucn_spp_range_area_%s_mol.csv', api_version)),
col_types = 'dd__') %>%
distinct()
taxa_ranges <- read_csv(file.path(dir_data,
sprintf('spp_marine_maps_%s.csv', api_version)),
col_types = 'ddcicccl') %>%
filter(comp_assessed) %>%
mutate(shp_base = basename(tolower(dbf_file)),
shp_base = str_replace_all(shp_base, '.dbf|_part.+', '')) %>%
left_join(spp_ranges, by = 'iucn_sid') %>%
left_join(shp_to_taxa, by = 'shp_base') %>%
filter(comp_assessed &!is.na(taxon)) %>%
mutate(log_range = log10(range_km2))
x <- ggplot(taxa_ranges, aes(x = log_range, ..scaled..)) +
ggtheme_plot() +
theme(panel.grid.major = element_blank(),
axis.text.y = element_blank(),
axis.title.y = element_blank(),
strip.text.y = element_text(angle = 0, hjust = 0)) +
geom_density(color = 'grey20', size = .25, fill = 'grey30') +
scale_x_continuous(limits = c(0, NA),
breaks = c(0:8),
labels = 10^(0:8)) +
labs(x = 'Range, km^2') +
facet_grid(taxon ~ .)
# ggsave(filename = file.path(dir_git, 'ms_figures', 'fig_SI_taxa_ranges.png'),
# width = 4.75, height = 3, units = 'in', dpi = 300)
print(x)